── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidylog)
Attaching package: 'tidylog'
The following objects are masked from 'package:dplyr':
add_count, add_tally, anti_join, count, distinct, distinct_all,
distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
full_join, group_by, group_by_all, group_by_at, group_by_if,
inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
relocate, rename, rename_all, rename_at, rename_if, rename_with,
right_join, sample_frac, sample_n, select, select_all, select_at,
select_if, semi_join, slice, slice_head, slice_max, slice_min,
slice_sample, slice_tail, summarise, summarise_all, summarise_at,
summarise_if, summarize, summarize_all, summarize_at, summarize_if,
tally, top_frac, top_n, transmute, transmute_all, transmute_at,
transmute_if, ungroup
The following objects are masked from 'package:tidyr':
drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
spread, uncount
The following object is masked from 'package:stats':
filter
library(RColorBrewer)library(viridis)
Loading required package: viridisLite
library(sdmTMB)library(sdmTMBextra)
Attaching package: 'sdmTMBextra'
The following objects are masked from 'package:sdmTMB':
add_barrier_mesh, dharma_residuals, extract_mcmc
library(patchwork)
Warning: package 'patchwork' was built under R version 4.3.3
library(RCurl)
Attaching package: 'RCurl'
The following object is masked from 'package:tidyr':
complete
library(tidylog)# devtools::install_github("seananderson/ggsidekick") # not on CRANlibrary(ggsidekick)theme_set(theme_sleek())# Set path:home <- here::here()
d <-read_csv(paste0(home, "/data/clean/temp_data_for_fitting.csv"))
Rows: 98616 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): area, source, db, id
dbl (8): year, temp, yday, month, day, VkDag, stn_nr, section_nr
date (1): date
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d <- d |>mutate(area =as.factor(area),source_f =as.factor(source),year_f =as.factor(year) ) |>filter(!area %in%c("VN", "TH")) # VN: no logger data, TH: to short time series
mutate: converted 'area' from character to factor (0 new NA)
new variable 'source_f' (factor) with 3 unique values and 0% NA
new variable 'year_f' (factor) with 83 unique values and 0% NA
filter: removed 5,545 rows (6%), 93,071 rows remaining
# Keep track of the years for which we have cohorts for matching with cohort datagdat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/clean/dat.csv")
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 10 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): area
dbl (1): mean_temp
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
order <- order |>mutate(area2 =ifelse(area %in%c("SI_HA", "BT"), paste0(area, "*"), area))
mutate: new variable 'area2' (character) with 10 unique values and 0% NA
group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 5,190 rows (2%), 333,270 rows remaining
filter (grouped): removed 107,058 rows (32%), 226,212 rows remaining
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
d <-left_join(d, gdat, by ="area") |>mutate(area =as.factor(area),growth_dat =ifelse(year >= min & year <= max, "Y", "N") )
left_join: added 2 columns (min, max)
> rows only in x 0
> rows only in y ( 2)
> matched rows 93,071
> ========
> rows total 93,071
mutate: converted 'area' from character to factor (0 new NA)
new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Drop data in SI_HA and BT before onset of warmingd <- d |>mutate(discard ="N",discard =ifelse(area =="BT"& year <=1980, "Y", discard),discard =ifelse(area =="SI_HA"& year <=1972, "Y", discard) ) |>filter(discard =="N")
mutate: new variable 'discard' (character) with 2 unique values and 0% NA
filter: removed 1,146 rows (1%), 91,925 rows remaining
# Drop heated areas actually for the full plotdf <- d |>filter(!area %in%c("BT", "SI_HA"))
nareas <-length(unique(order$area)) +2# to skip the brightest colors that are hard to readcolors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)[-c(6, 7)]
Plot data
df |>filter(growth_dat =="Y") |>distinct(area, year, source) |>ggplot(aes(year, area, color = source)) +geom_point(position =position_dodge(width =0.4), shape =15) +labs(y ="Site", x ="Year", color ="Source") +scale_color_brewer(palette ="Dark2") +theme(legend.position ="bottom")
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
ggsave(paste0(home, "/figures/supp/qq_temp_area.pdf"), width =17, height =17, units ="cm", device = cairo_pdf)# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclearnd_area_sub <- nd_area |>mutate(keep ="N",keep =ifelse(area =="FM"& year <=1980, "Y", keep), # use FM instead of BTkeep =ifelse(area =="SI_EK"& year <=1972, "Y", keep) ) |># use SI_EK instead of SI_HAfilter(keep =="Y") |># Now change the labels to BT and SI_EK...mutate(area =ifelse(area =="FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modellingnd_area <-bind_rows(nd_area, nd_area_sub) |>select(-keep) |>mutate(growth_dat =ifelse(area =="SI_HA"& year %in%c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)
# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by yearnd_area |>mutate(area =ifelse(area %in%c("SI_HA", "BT"), paste0(area, "*"), area)) |>ggplot(aes(yday, y = year, fill = pred, color = pred)) +geom_raster() +facet_wrap(~area, ncol =4) +scale_fill_viridis(option ="magma") +scale_color_viridis(option ="magma") +labs(x ="Day-of-the-year", y ="Year", color ="Predicted SST (°C)", fill ="Predicted SST (°C)") +theme(legend.position =c(0.78, 0.14))
mutate: changed 182,268 values (20%) of 'area' (0 new NA)
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
ggsave(paste0(home, "/figures/supp/temp_pred_yday_area.pdf"), width =17, height =17, units ="cm")ggsave(paste0(home, "/figures/for-talks/temp_pred_yday_area.pdf"), width =14, height =14, units ="cm")for (i inunique(nd_area$area)) { plotdat <- nd_area |>filter(area == i)print(ggplot(plotdat, aes(yday, pred, color = source)) +scale_color_brewer(palette ="Dark2") +facet_wrap(~year) +geom_point(data =filter(d, area == i & year >min(plotdat$year)), size =0.2,aes(yday, temp, color = source) ) +geom_line(linewidth =0.3) +labs(title =paste("Site = ", i), color ="", linetype ="") +guides(color =guide_legend(title.position ="top", title.hjust =0.5)) +theme_sleek(base_size =8) +theme(legend.position ="bottom",legend.direction ="horizontal",legend.spacing.y =unit(-0.3, "cm") ) +labs(x ="Day of the year", y ="Predicted SST (°C)") )ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_", i, ".pdf"), width =17, height =17, units ="cm")}
mutate_at: converted 'lat' from character to double (0 new NA)
converted 'lon' from character to double (0 new NA)
order <-read_csv(paste0(home, "/output/ranked_temps.csv")) |>mutate(area2 =ifelse(area %in%c("SI_HA", "BT"), paste0(area, "*"), area))
Rows: 10 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): area
dbl (1): mean_temp
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: new variable 'area2' (character) with 10 unique values and 0% NA
preds_area |>mutate(area2 =ifelse(area %in%c("SI_HA", "BT"), paste0(area, "*"), area)) |>ggplot(aes(year, temp, color = temp)) +facet_wrap(~factor(area2, levels = order$area2), ncol =5) +geom_line(linewidth =0.6) +labs(x ="Year", y ="Model-predicted annual average temperature (°C)") +scale_color_viridis(option ="magma", name ="Site") +guides(color ="none")
mutate (grouped): new variable 'area2' (character) with 10 unique values and 0% NA
ggsave(paste0(home, "/figures/annual_average_temperature.pdf"), width =17, height =9, units ="cm")# Check year rangepreds_area |>group_by(area) |>summarise(min =min(year),max =max(year) ) |>arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
area min max
<chr> <dbl> <dbl>
1 JM 1953 2016
2 BT 1963 2000
3 FM 1963 2000
4 SI_HA 1966 2014
5 SI_EK 1968 2015
6 FB 1969 2016
7 MU 1981 2000
8 HO 1982 2016
9 BS 1985 2002
10 RA 1985 2006
# Check temperature rangepreds_area |>group_by(area) |>summarise(min =min(temp),max =max(temp) ) |>arrange(min)
group_by: one grouping variable (area)
summarise: now 10 rows and 3 columns, ungrouped
# A tibble: 10 × 3
area min max
<chr> <dbl> <dbl>
1 RA 3.06 9.56
2 JM 3.37 11.8
3 HO 3.99 8.07
4 FB 5.04 10.6
5 MU 5.16 8.57
6 SI_EK 5.49 10.4
7 BT 5.87 16.6
8 FM 5.87 15.0
9 BS 6.22 11.8
10 SI_HA 7.76 16.9
# Save prediction dfwrite_csv(preds_area, paste0(home, "/output/gam_predicted_temps.csv"))
